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Abstract— Some recent results obtained using solution-adaptive finite element methods in two-dimensional 
problems in linear elastic fracture mechanics are presented. The focus is on the basic issue of adaptive 
finite element methods for validating the new methodology by computing demonstration problems and 
comparing the stress intensity factors to analytical results. 


1. INTRODUCTION 

One of the most difficult analytical challenges in 
engineering mechanics is the modeling of flawed 
structures and the computation of the structural 
response and flaw propagation. These models are 
generally quite large for many space vehicle struc- 
tures. The commercially available codes are often 
used to model the unflawed structure and to perform 
a stress analysis. These large models must then be 
reconstructed manually to introduce a flaw and to 
remesh the problem for accurate modeling of the 
flaw. One of the most technically challenging areas is 
fatigue crack propagation in these structural com- 
ponents. Specifically, fatigue life prediction error 
sources such as the stress intensity factors (SIF) are 
difficult to quantify. Practically, SIFs are calculated 
either from handbooks or other simplified equations. 
However, these equations are applicable only for 
structural components that approximate the defined 
numerical and geometric conditions. Furthermore, 
handbooks do not usually consider mixed mode 
crack situations which are common for actual struc- 
tures. Thus, there is a need to minimize the error 
associated with SIF values and have them integrated 
into fatigue models for improved fatigue life 
predictions. 

Much of contemporary fracture analysis still fol- 
lows classical ideas. Modeling of complex geometries 
has been difficult and conventional computational 
procedures cannot provide accurate stress predictions 
for bodies of complex shape. Many fracture theories 
are developed with only simple stress states in mind 
and are seldom factored into realistic stress environ- 
ments of the type actually experienced in working 
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structural components. In concept, many of the 
shortcomings of classical fracture analysis may be 
overcome through the development of more sophisti- 
cated computational models. By using new finite 
element capabilities and new concepts in fracture 
mechanics, it should be possible to study a variety of 
basic issues connected with fracture and crack 
growth. These include the use of more elaborate 
models of material constitution and component 
geometry, more general crack initiation and growth 
criteria, more accurate methods for prediction crack 
development, and more physically reasonable models 
of crack arrest mechanisms. Therefore, an effort has 
been made on the development of the efficient and 
accurate modeling technique of large space structures 
containing flaws using solution-adaptive methods. 
These procedures, using finite element methods, auto- 
matically adjust the grid points for refinements of 
meshes of quadrilateral elements to produce a mini- 
mum error solution. A viable approach is to develop 
a new fracture mechanics analysis tool which is based 
on modern adaptive finite element methods. The 
primary goals of this study are to: (1) develop an 
advanced and reliable numerical method for perform- 
ing linear elastic analysis of flawed structural com- 
ponents; (2) validate the new methodology by 
computing demonstration problems and comparing 
the stress intensity factor to analytical results. 

As a result, efficient solution-adaptive algorithms 
are derived which are suitable for large-scale compu- 
tations for the solution of confined crack regions and 
a two-dimensional fracture mechanics analysis code 
is developed in which a discrete least-squares algor- 
ithm and an energy release method are implemented 
and validated. Demonstration problems are also 
computed and the SIFs compared with analytical 
formula. The agreement is reasonably good and thus 
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provides validation of the new methodology. A 
general procedure and methods employed are given 
in Secs 2 and 3, and the results of demonstration 
problems are given in Sec. 4. 

2. ADAPTIVE COMPUTATIONAL METHODS 
2.1. Adaptive schemes 

Adaptive methods, which are numerical schemes 
which automatically adjust themselves to improve 
numerical solutions, are generally composed of 
three key ingredients. The first of these is an 
adaptive scheme which is capable of dynamically 
modifying the structure of the computational grid. 
The second crucial component is an error estimator 
or error indicator which is used to identify regions 
of the computational domain which contain 
relatively large numerical errors. The final necessary 
ingredient is an adaptive methodology which 
combines the mesh modification schemes with the 
error estimators to produce a robust and efficient 
refinement sequence. 

The first work on adaptive finite element methods 
was presented in 1971 by Oliveira [1] which discussed 
grid optimization by minimizing the energy through 
optimal node redistribution. This type of approach, 
node redistribution, is the basis of moving mesh 
adaptive methods (r -methods) developed for both 
solid mechanics problems and flow analysis [1-8]. 
Following this work, several other adaptive finite 
element schemes have emerged which include h - 
refinement methods, which subdivide the elements 
of the computational domain [9-13], p -enrichment 
methods, which increase the polynomial order of 
the approximations [14-18], and h-p methods, 
which increase the polynomial order, as well as 
subdivide the elements comprising the grid [19-23]. 
In these schemes, the mesh is automatically refined 
when the local error indicator exceeds a preassigned 
tolerance. The h -scheme presents a difficult data 
management problem, since they involve a 
dynamic regeneration of the mesh, renumbering 
of grid points, cells or elements, and element 
connectivities as the mesh is refined. Since it has 
been proven that the h -methods can be very 
effective in producing near-optimal meshes for 
given error tolerances and also be used to coarsen 
a mesh (use large mesh cells and thereby reduce 
the number of unknowns) when the local-error 
becomes lower than an assigned lower-bound 
tolerance among the adaptive methods [9], the h - 
method is taken in this study. A sample calculation 
obtained with our h -method is shown in Fig. 1 for 
a uniformly loaded plate containing a crack. Our 
procedure dynamically refines the mesh, assigning 
large elements where the error is small and small 
elements where the error tends to be large, thus 
capturing the singularity of the crack tip. Computed 
contours of the a y stress component are also shown 
in the figure. 




Fig. I. Uniformly loaded crack plates: h -adapted mesh and 
o y stress contours. Note that the model is not exactly 
symmetric for the top and bottom parts. 


2.2. Error estimation 

An adaptive scheme is essentially useless without a 
reliable, accurate, and efficient error indicator. The 
error indicator provides an estimate of the local error 
in the numerical solution and an indication where the 
grid should be adapted. In general, there are two 
distinct types of error estimators, interpolation error 
estimators and residual error estimators. The inter- 
polation error estimates typically provide only a very 
crude estimate of the error, but usually provide a very 
good indication of the relative errors and are quite 
inexpensive to implement. Residual methods, on the 
other hand, typically provide a much better estimate 
of the actual error in the numerical solution but are 
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generally much more computationally expensive. 
Therefore, the interpolation error estimates are taken 
in this study. Some general properties of the distri- 
bution of error in finite element meshes and the use 
of so-called interpolation error estimates shall be 
discussed here. This easily implementable local esti- 
mator sometimes provides only a rather crude 
estimate of the actual local error but can be divised 
to give a correct indication of the relative error 
between successive meshes or approximation orders 
and, thus, correctly direct an adaptive strategy to 
systematically reduce local error. For a more detailed 
discussion of interpolation theory, see [24, 25]. The 
simplified theory of this error estimates follows. 

Let u be a smooth function defined over a regular 
domain fic 5R 2 . The lF rp (ft)-seminorm of u is defined 
by 





d ,+J u p 
3x ] dx{ 


U* 0 



( 1 ) 


Suppose that u on the right-hand side is now replaced 
by a finite element approximation u h and 
kLk + i., = l*4 + i , P + 0(h). Then the above estimate 
indicates that the local error in the W m ' p (G )-semi- 
norm is proportional to the error indicator, 
^ («//>' )-(«/*>) + * + i - | A t + , Some choices for n, m, k, p , 

and p' of interest are: 

(i) n = 2, m = 0, k = 1, p = p' = 2 then 

\E h \ LHG) <C -h 2 -\u\ 1XG . ( 7 ) 

In this case, one must approximate the W /2,2 -semi- 
norm of u over ft; i.e., the L^-norm of the second 
partial derivatives of u. The error indicator 0 e is then 
set equal to \E h \ L2{(ie) for finite element ft,. 

(ii) n =2, p = oo, p' = 1, k =0, m =0 then 

| E h \ LHG) = C ' h 2 • l^veragel ^ C ' h* ' |w|i, x .. G 

^ C ■ A 3 ■ max |V • w(x)| (8) 

.veG 


where 1 < p < go and r is a non-negative integer. For 
the special case ofp = oo, which is also of interest, the 
IF'^ftJ-seminorm is given by 


“ max max 

ij > 0 

With these definitions of the seminorm, the Sobolev 
norm of u is then 

II M II W'-Hfl) = ^ X 1^1 ^ 

Let G be an arbitrary convex subdomain (a finite 
element) of ft over which u is interpolated by a 
function u h which contains complete piece poly- 
nomials of degree k. Then it can be shown that 
the local interpolation error in the H /rM,p (n)-semi- 
norm is 

^ c ■ — ■ h (n/p) - • |«l»*+' P( 0 ). (4) 

y m 


d i+J u(x) 


dx\ dx{ 


( 2 ) 


then we have for the error indicator 0 e 

\E h L„ lge „v„ Q ,*e e = h- max |V • u (x )| . (9) 

In all of these cases, it is also possible to estimate the 
constant C. While we shall not describe how this is 
done (see [25]), our experience is that it is a worth- 
while computation that can lead to schemes with 
good effectivity indices. 

2.3. Adaptive methodology 

We shall now suppose that an error indicator 6 e can 
be calculated for each finite element Q e in a given 
mesh. The error indicator is, in general, a real number 
representing the local error on a suitable norm, and 
it is computed using one of the procedures described 
in the next section. 

The decision to refine the mesh is based on whether 
or not local-indicators exceed preassigned tolerances 
and can be summarized by the following steps. 

h -Refinement / unrefinement methodology. The h- 
procedure involves the following steps: 


where h is the diameter of domain G, y is the diameter 
of the largest sphere that can be inscribed inside G, 
n is the dimension of the domain ft, C is a constant 
independent of h, and p' is pi(p — 1)- If 7 is P r0 ~ 
portional to h , and if it remains proportional in 
refinements of G defined by parametrically reducing 
A, we have 

\E h L^c < c * h (nip,) ' {n!p)+k + i - m ■ + (5) 

with 

I ' \m,p\G = I * \w^P\G)y etC> aiK * ^ =U—U h .( 6 ) 


1. For a given domain ft, a coarse finite element is 
constructed which contains only a number of ele- 
ments sufficient to model basic geometric features of 
the domain. 

2. As the adaptive process is designed to handle 
groups of four elements at a time, a finer starting grid 
may be generated by a bisection process, if desired, 
to obtain an initial set of element groups. 

3. The numerical solution is calculated on this 
initial coarse grid, and the error indicators 6 e are 
computed over all M elements in the grid. Let 

0 max = max 0 t . (10) 
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4. The groups of elements are scanned and the 
group errors are computed 

Coup = X 0,*, (11) 

k = I 

where e k is an element number in group k and p — 4. 

5. Error tolerances are defined by two real num- 
bers, 0 < ot, < 1. If 

0max (12) 

element e is refined. This is done by bisecting element 
e into four new subelements. If 

Coup « a ' C, (13) 

group £ is unrefined by replacing this group with a 
single new element with nodes coincident with the 
corner nodes of the group. 

This general process can be followed by any choice 
of an error indicator. Moreover, it can also be 
implemented with any boundary frequency in the 
solution of transient problems or problems with 
time-dependent boundary conditions. 

3. CALCULATIONS OF SIFs 

Fracture mechanics within the linear theory of 
elasticity is known as linear elastic fracture mechanics 
(LEFM). Linear elastic fracture analysis is valid if the 
assumption of small-scale yielding, namely when the 
inelastic region in the vicinity of each crack tip is well 
within the range of validity of the asymptotic sol- 
utions of the crack in an infinite elastic medium. 
Assuming that the crack faces are free of traction, the 
asymptotic expressions for the neighborhood of the 
crack tip are given in [34]. In this study the stress 
intensity factors were calculated using one of the 
following methods. 

3.1. Discrete least -square fits 

Consider a patch of elements around the crack tip 
as shown in Fig. 2; for example, we may define 

ftiip = {* e QJ \x™' - x l * p | ^R(h)}. (14) 

Let us define the discrete least squares functional 

nintp 

J(K ly K u )= X (<r l-eyny 

j - 1 

X D‘ x '(<!*; (15) 

Here, xf nt is the position vector for the centroid of 
element Q e , x tip is the position vector for the tip of 
the crack, R(h) is the radius of the disk which 
includes the centroids of the element of the subdo- 
main, which is assumed to be a function of the mesh 
size /t, a 1 } ~ e h (Xj) denotes the approximate stress 
vector at point x Jf a* 55 ™ = Aj, K u ) 



Fig. 2. A subdomain which includes the elements of the 
grids with centroids in the interior of the disk 
{«* e — „v| /?(/?)}. The elements which include the 

crack tip appear shaded. 

denotes the asymptotic stress vector at point Xj, D~ l 
is the inverse of the material modulus matrix, x } 
denotes an investigation point from the set of inte- 
gration points of the elements which form ft up , W J 
denotes the weight for the j th integration point, and 
nintp denotes the total number of integration points 
in the subdomain Q tip . 

The components of the stress a asym (x, K x , A n ) de- 
pend linearly on the SIFs; their analytical expression 
is given in [33]. Approximate values for the SIFs are 
extracted as minimizers of the discrete least squares 
functional J{K { , K l{ ); the condition of a stationary 
value leads to a two-by-two system of linear 

equations for the approximate values of K { and K n , 

namely 

Cj ^ °mtp ^(^sym)r 

dK x ~ ~ dK { 

x D~ l • (<Tj — <Tj sym * W } ) = 0 (16) 

dJ ^ n ^p^((r; sym ) r 

dK n ~~ dK n 

xD-'-(<r*-eW m 'Wj) = 0. (17) 

A bilinear element is employed to compute the ap- 
proximate finite element solution and the discrete 
least squares functional in each element is defined 
using 1 x 1 and 2x2 Gauss integration rules. Let K\ , 
and K \ , and K\ ,, denote the extracted values of 
the SIFs corresponding to the use of the 1 x 1 and 
2x2 rules in the definition of the functional. Since 
the definition of the K 1 values involves superconver- 
gent values of stress at the centroids of the elements, 
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(perpendicular to 



Fig. 3. Planar crack in three dimensions. 


the following hypothesis was taken from the numeri- 
cal experiments to determine that the error in the K 
values of the SIFs can be estimated by comparing 
them with the K ] values 

K-K 2 ^K'-K 2 . (18) 

3.2. Equivalent domain integral methods 

Equivalent domain methods were developed by 
Shih, Moran and co-workers [26-28]. Li et ai [26] 
derived a volume integral expression for the energy 
release rate which can be naturally employed to 
extract the energy release rate along three- 
dimensional crack fronts from three-dimensional 
finite element calculations. Also, the two-dimensional 
area integral analog of the volume integral expression 
was described. The energy release rate is defined 
as follows: Consider a planar crack front with a 


continuously turning tangent as shown in Fig. 3. Let 
<5 L(j) denote the crack advance at the points in the 
direction normal to the crack front and d5 denote the 
length of a crack element along the crack front. Then, 
within first-order terms in <5L(s) and (7(5), namely the 
pointwise energy release rate at the location s of the 
crack front, defined by (see [26] and references 
therein) 

f (7(5) ’ SL(s) ds = — Sn. (19) 

J crack front 

Here Sn denotes the change in the total potential 
energy due to the local advance of the crack given by 
SL(s). The above definition may be considered as a 
weak or variational statement satisfied by the energy 
release rate function, which is defined along the crack 
front. Li et ai [26] suggested a volume-integral 


% (x) 
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Fig. 5. Definition of the parameters involved in the definition of the components of the 7-integral. 


expression for the energy release rate; the two- 
dimensional version of this expression reads as 
follows (see [26], pp. 408^409 for the derivation): 


G =J = 





( 20 ) 


Here A l denotes the simply connected region of 
integration in Fig. 6 and 


du ( 

is the Xj-component of Eshelby's energy momentum 
tensor [29, 30], W is the strain energy density, and q { 
is a sufficiently smooth function in A { which is unity 
on the inner boundary C ] and vanishes on the outer 
boundary C 2 of A { . 

A general study of crack-tip contours and domain 
integrals was presented by Moran and Shih [27, 28]. 
In two dimensions, the energy release rate is ex- 
pressed by an equivalent domain integral given by 
eqn (20); a schematic definition of the integration 
domain A { and the weight function q l is given in 
Fig. 4. According to [31] the basic definition of 



Fig. 6. Simply connected region A { with subsets C, and C 2 
of the boundary. 


^-integral components, as crack-tip parameters, is 
(see Fig. 5) given by 

(2i) 

Here W is the stress-work density, a t) are stresses, u, 
are displacements, n k are components of the unit 
normal vector at points on the contour f. A local 
crack front coordinate system is given by x, and x 2 : 
x, is normal to the crack front, x 2 is orthogonal to . 
In the linear elastic case, it is possible to define SIFs 
through the calculation of J k (k = 1,2). According 
to [32] the relationships between the /-integral com- 
ponents and the stress factors are 


* I = - Jl + y/Ji+Ji) (22) 

*11 = - yfj\ + / 2 ). (23) 

where 

f E 

E* — J | v 2 f° r plane strain state 

for plane stress state. 
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Fig. 8. Edge cracked beam. 


4. NUMERICAL examples 4.1. Single edge crack ( problem 1) 

Figure 7 is a schematic of this demonstration 
This section presents the demonstration problems. problem The problem parameters are given as 

The cases consist of a single edge crack and an edge f 0 n 0ws Dimensions: L = 12 inches, W = 4 inches, 

cracked beam; both cases were computed using two and T= , oinch. Material properties: E = 1 0 6 psi 

different crack lengths. The results for each case are and y _ 0 3 Crack lengths: a = 2.0 and 0.2 inches, 

shown as stress contour plots on the deformed shape Stress d i str ibution: a = 100 ksi (top and bottom), 

scaled by a factor of 100 to clearly show crack. Mode I stress intensity factors were calculated at each 

Fifteen to twenty levels of refinement were used to crack i engt h f or the given stress distributions and 

compute the SIFs using the methods shown in Sec. 3. both crack | engt h s f or two cases Q f this problem. 

The ‘equivalent-domain integral' method was used 

for problems 1 and the ‘discrete least-squares fits’ 4.1.1. Crack length a = 2.0 (problem \A ). The rep- 
method was used for problem 2. resentative solutions for this problem are shown in 



Fig. 9. Problem 1A: (a) initial grid, (b) fourth level refinement grid. Note that the grids are scaled by a factor of 100 to 

clearly show the crack. 
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Table 1. SIF comparisons for example problem 1 



Adaptive FEM 

Ref. [33] 

Ref. [34] 

Problem 1A 

715 

755 

706 

Problem IB 

87.5 

89.9 

89.6 


Figs 9 and 10. Figure 9 shows the initial grid distri- 
bution and fourth level refinement. Figure 9 shows 
the grid with a x and a y contours for the fourth level 
of refinement. 

4.1.2. Crack length a = 0.2 ( problem 15). The rep- 
resentative solutions for this problem are shown in 
Figs 11 and 12. Figure 11 shows the initial grid 
distribution for this case. Figure 12 shows the grid 
and fourth level refinement with a x and a y contours 
for the fourth level of refinement. 

4.1.3. Validation for problem 1. The ‘energy 
release’ method was used in this problem. 
References [33] and [34] provide formulas for 
estimating the SIF for this problem. The formulas 
used are of the form 

K x = a ■ \J 7t • a ■ F(h), h = (24) 

The references use fourth-order expansion in h to 
determine the multiplier F. Table 1 compares K x 
obtained from the adaptive finite element code and 
the analytical formulas from the references. 



4.2. Edge cracked beam ( problem 2) 

Figure 8 is a schematic of this problem. The 
problem parameters are given as follows. 
Dimensions: W = 2 inches, L — 41V = 8 inches, T = 
1.0 inch. Material properties: £ = 10 6 psi, v=0.3. 
Crack lengths: a =0.5 and 0.1 inches. Stress distri- 
bution: P = 5000 lb (point load). Mode I stress 
intensity factors were calculated for both crack 
lengths for two cases of this problem. 

4.2. 1 . Crack length a = 0.5 ( problem 2 A ). Figure 1 3 
shows the initial and fourth level refinement grid 
distribution. The solutions of fourth level of refine- 
ment for this problem are shown in Fig. 14 to show 
the grid with o x and a y contours for this level of 
refinement. 

4.2.2. Crack length a = 0.1 ( problem 25). The sol- 
utions of this problem are shown in Figs 1 5 and 1 6. 
Figure 16 shows the grid with a x and a y contours for 
the fourth level of refinement. 

4.2.3. Validation for problem 2. The ‘least squares’ 
method was used in this problem. References [33] and 
[34] provide formulas for estimating the SIF for this 
problem. The formula used are of the form 

K y = a jn - a ■ F(h), h=^- 



Fig. 12. Problem IB: (a) a, for fourth level refinement grid, (b) a for fourth level refinement grid. Note that the grids are 

scaled by a factor of 100 to clearly show the crack. 
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Fig. 13. (a) Initial grid; (b) fourth level refinement grid of problem 2A. 





Fig. 14. (a) a x on the fourth level refinement grid; (b) <r, on the fourth level refinement grid of problem 2A. (Note: grid 

scaled by a factor of 1 00 to clearly show the crack.) 



443 


I 



of problem 2B. 
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Table 2. SIF comparisons for example problem 2 



Adapative FEM 

Ref. [33] 

Ref. [34] 

Problem 2A 

17,901 

18,917 

20,301 

Problem 2B 

8611 

8418 

8985 


The references uses the fourth-order expansions in h 
to determine the multiplier F. Table 2 compares K t 
obtained from this study with K x obtained using the 
analytical formula from the references. 

5. CONCLUSIONS 

A solution-adaptive computational method has 
been developed and implemented for accurately de- 
termining the SIFs of two-dimensional crack 
configurations with predetermined crack locations 
and loading conditions. This new method and pro- 
cedures have been validated for demonstration prob- 
lems showing reasonably good agreement with 
analytical formula and can be used for performing 
the two-dimensional fracture mechanics analysis in 
complex structures. Also, this study demonstrated the 
feasibility of the methodology developed to be ex- 
tended in the three-dimensional analyses. 
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